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We reexamine whether the essence of high-T c superconductivity is contained in doped Hub- 
bard models on the square lattice by using recently developed pre-projected Gaussian-basis 
Monte Carlo method. The superconducting correlations of the d x i_ y i wave symmetry in the 
ground state at distance r decays essentially as r -3 . The upper bound of the correlation at 
long distances estimated by this unbiased method is 10 -3 , indicating that recent extensions 
of dynamical mean-field theories and variational methods yielded at least an order of magni- 
tude overestimates of it. The correlations are too weak for the realistic account of the cuprate 
high-T c superconductivity. 
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Since the proposal that the simple Hubbard model on 
a square lattice captures essential physics of the high- 
temperature copper-oxide superconductors, 1 ) the model 
has been extensively studied and now become one of 
the most intensively studied issues in condensed mat- 
ter physics. However, the ground state of this model, 
particularly whether the high-T c superconductivity is ac- 
counted for, has long been highly controvert ial. 

The Hubbard Hamiltonian contains the on-site 
Coulomb interaction U and the hopping integral t as 



n s 



n s 



H 



(1) 



where c\ a (ci a ) is the creation (annihilation) operator of 
an electron with spin a on the z-th site for 7V s -site lattice 
and (i, j) represents nearest neighbor pairs. 

In the weak to intermediate coupling region, the 
auxiliary-field quantum Monte Carlo (AFQMC) method 
gave negative results for the superconductivity around 
up to U/t = 4 by the analysis of the superconducting 
correlation functions. 2, 3 ) Zhang et al. calculated pairing 
correlations of the 4 2 -y 2_wave channels by combining 
the constrained-path approximation (CPMC) and con- 
cluded the absence of the superconductivity. 4 ) Due to 
the negative sign problem known as a major obstacle of 
the Fermion simulation methods, the AFQMC method 
does not offer converged results for larger U/t. 

On the other hand, with the dynamical cluster ap- 
proximation (DC A), Maier et al. claimed, from the clus- 
ter sizes up to 26 sites, that the d £C 2_ 2/ 2-wave pair-field 
susceptibility indicates the transition to the d x 2_ y 2-wdive 
superconductor at the critical temperature T c « 0.023t 
at 10% doping, even at U/t = 4. 5 ) Cellular dynamical 
mean- field theory (CDMFT), cluster perturbation the- 
ory and variational cluster theory (VCT) gave similar 
results. 6 8 ) Accuracy of the mean- field treatment con- 
tained in these studies remains to be examined by larger 
clusters. Numerical results based on variational Monte 



Carlo (VMC) methods suggest the firm d x 2_ y 2-wdive su- 
perconductivity for U/t > 6 upon doping. 9 ) However, 
the VMC method could overestimate the stability of su- 
perconductivity by relatively worse estimate of energy 
for correlated metals. From these controversies, one finds 
that the stability of superconductivity in Hubbard mod- 
els still remains a fundamental open issue. In the cir- 
cumstance of the ground state highly competing with 
others, unbiased accurate approaches are needed to an- 
swer whether the realistic amplitude of the copper-oxide 
superconductivity is accounted for. 

In this letter, we reanalyze the pairing correlations of 
the d x 2_ y 2 symmetry on doped Hubbard models by using 
recently developed pre-projected Gaussian-basis Monte 
Carlo (PR-GBMC) method. 10 ) Thanks to the absence of 
the negaitive sign problem, this method allows us to go 
beyond the tractable range of the conventional methods 
including doped and larger U region without any approx- 
imation. We reveal that the pairing correlations do not 
show distinct enhancement up to U/t = 7. An accurate 
estimate of the upper bound of the superconducting cor- 
relation (~ 10 -3 ) is far below the recent results by DCA, 
CDMFT, VMC and VCT, indicating the necessity of the 
higher accuracy than these studies in the literature for 
assessing the possible superconductivity in the Hubbard 
models. 

Our PR-GBMC method is described in detail and ex- 
tensively benchmarked elsewhere. 10 ) Here we only ex- 
plain the basic procedure of PR-GBMC. To calculate the 
ground state of the Hamiltonian (1), we solve the Liou- 
ville equation of the density-matrix operator p{r) (here, 
r denotes the imaginary time): 

dp(r) _ 1 
6V 2 

using an expansion of p(r) by a Gaussian-basis A(fi, n): 
p( T ) = j dfldnP(fl,n;T)A(fl,n), (3) 



H,p(r) 



(2) 
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A(^,n) = ^det(/-n):e-" t [ 2/+ ^ T - / ) ^ :, (4) 

where P(fi, n; r) is expansion coefficient, J is a 2N S x 2N S 
unit matrix and the vector operator c) is defined by 



(cjpcjp--- , c]y s |, c{|, • • • ,c]y s |). (5) 

The colon bracket : : is the normal ordering oper- 
ator. One of the parameters of the Gaussian-basis CI 
works as the weight of the importance sampling of 
the Monte Carlo procedure and the other parameter 
n is a 27V S x 2N S matrix characterized by ^(ia),(ja) = 
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Tr 
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/CI. By using the ex- 



pansion (3), the Liouville equation (2) is mapped into 
a Fokker-Planck equation for P. It should be noted 
that P can be taken positive definite. 10 ' 11 ) In the ac- 
tual calculation, we solve the corresponding Langevin 
equations with respect to CL and n. 12 ) In our PR-GBMC 
method, 10 ^ by using quantum- number projectors 13 ^ P = 
/ dxg(x)e i& h ^ d , the weight CL is transformed into a 
quantum- number-projected weight Cl. Here, h(x) is a 
2N S x 2N S Hermitian matrix and g(x) is an integration 
weight with x being a parameter of the quantum- number 
projection such as phase, Euler angles of the spin-space, 
etc., . The importance sampling is performed by the 
quantum-number-projected weight CL defined by 



Cl = CL dxg(x) det (e 



Jh(x) 



J)n 2 



(6) 



We have taken large enough r to assure the convergence 
to the ground state and confirmed the absence of the 
boundary term in P following the procedure in Ref.. 10 ) 

To reduce finite-size effects and estimate it in a sys- 
tematic fashion, we mainly employ fillings at which the 
corresponding noninteracting systems show closed-shell 
structure. 3 ^ For comparison, open-shell conditions are 
also studied in some cases. 

The accuracy of the present method is seen typically 
in the estimate of the energy per site for the Hubbard 
model at half filling: 10 ) -0.8575 ± 0.0005 for 6 x 6 lattice 
in comparison with —0.8575 ± 0.0008 by the accurate 
AFQMC result and -0.8595 ± 0.0005 for 8 x 8 lattice in 
comparison with -0.8607±0.0009 by the AFQMC result. 
This accuracy is more than one order of magnitude better 
than typical VMC accuracy. 

We now present results for the equal-time pairing cor- 
relation defined as 



1 



N 



calculate the superconducting correlation defined by the 
summation of P a (r) over r: 



P «W = 2^ 2J A ^)A a (* + r) + A a (i)Al(i + r)>, (7) 
where A a is the order parameter defined as 



A Q (i) = i=^/ Q (r)(c j 



^ c i+rl 



Here, / a (r) is the form factor defined as 



(8) 



f2d(r) = 5r y ,o(3r x ,l + ^-l) - 5r x ,o(5r y ,l + *r v ,-l) (9) 

with Sij being Cronecker's delta. The suffix a — 2d rep- 
resents the cl^^ 2 -wave. In addition to P a (r), we also 



(10) 



where we multiply the factor 1 /4 to allow direct compar- 
ison with the data in Refs. 3 and 14. 
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Fig. 1. (Color online) Pairing correlation for d x 2 _ y i -wave as 
function of on-site interaction U/t. Squares and triangles rep- 
resent 6x6 lattice at the doping concentration S = 0.28 and 
8x8 lattice at 5 = 0.22, respectively. 



Figure 1 shows that S2d decreases first, then increases 
above U/t ~ 5 with the increase of U. The sharp en- 
hancement of S2d seems to be consistent with the en- 
hancement of the condensation energy above U/t ~ 6 
seen in the VMC results (see Fig. 3 in Ref. 9). The most 
of the enhancement observed in S^, however, comes from 
the short-range correlations. Figure 2 shows the pairing 
correlation of the d a .2_ 2/ 2-wave at U/t = 6 subtracted by 
that at U/t = 0. It shows that the on-site and nearest- 
neighbor pair correlations exhaust most of the enhance- 
ment in S2d- The enhancement in the short-range part 
was also pointed out by the CPMC results. 4 ) 

The enhancement of the on-site correlation P2d(r = 0) 
comes from the nearest-neighbor spin and charge cor- 
relations. In fact, P2d( r — 0) contains the sum over 
the nearest-neighbor spin and charge correlations, C = 
(—ASo-Sx+NqNz), where x denote the unit vectors in the 
x directions. S and N are spin and charge operators, re- 
spectively. The enhancement is well accounted for by the 
enhancement of C, which is not related to the supercon- 
ductivity but to enhanced antiferromagnetic correlations. 

Next, we analyze the filling dependence of the d x 2_ y 2- 
wave pairing correlation p2d(r). As we see in Fig. 3, the 
pairing correlation P2d(r) has the largest values at the 
doping S = 0.22. Although S = 0.22 satisfies the closed 
shell condition, the other two fillings are under the open 
shell condition. We do not find an appreciable difference 
between open and closed shell conditions. 

To gain insight into the main issue whether the non- 
zero offset exists in the long-range part of the pairing 
correlation, we analyze system size dependence of P2d(r) 
at the filling around 5 ~ 0.2, where the long-range part 
of P2d(V) may have the largest values. Figure 4 shows 
P2d{r) for several system sizes, where P2d(V) decreases 
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Fig. 2. (Color online) Spatial dependence of pairing correlations 
of the d x 2_ y 2-wave at U/t = 6 subtracted by those at U/t = 0. 
The triangles, circles and the squares represent sizes 6x6 lattice 
at the doping concentration 5 = 0.28, 8 x 8 lattice at 5 = 0.22 
and 10 x 10 lattice at 5 = 0.18, respectively. 
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Fig. 3. (Color online) Filling dependence of pairing correlation 
P2d(r) as function of distance r on the 8x8 lattice at U/t = 6. 
Squares, triangles and circles represent P2d( r ) at 5 = 0.34, 0.22 
and 0.09, respectively. 



with the increasing system size as well as with the in- 
creasing distance r at large r. It is known that in the non- 
interacting case, the pairing correlation decays asymptot- 
ically as P2d( r ) °c r ~ 3 as is inferred from 20x20 lattice 
result in Fig.4. 3 ) This behavior is gradually reproduced 
with the increase of the system size even at U/t = 6. It 
should be noted that in addition to the power-law de- 
cays, the correlations show increase again at the farthest 
distances in finite size systems (e.g., r > 3 for 6 x 6, 
r > 4 for 8 x 8 and r > 5 for 10 x 10 lattices in Fig. 4). 
These increases are well accounted for by the superim- 
pose of the tail from the other directions due to the pe- 
riodic boundary condition. After correcting this artifact, 
the pairing correlation shows no signal of level off to a 
nonzero correlation at long distances. From converged 
values of correlations up to r ~ 5, we may safely con- 
clude that \P2d(r — > oo)| < 10 -3 for the Hubbard model 
at U/t = 6. This main result of this paper poses a severe 
constraint on the possible superconductivity. 

When we could assume the existence of the d x i_ y i- 
wave superconducting long-range order, the order param- 
eter average (A2d) should be related to P2d as \P2d(r —> 
oo)| ~ (A 2 d) 2 . Although recent DCA, 5 ) CDMFT, 7 ) 
VCT 8 ) and VMC 9 ) studies are not necessarily consistent 
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Fig. 4. (Color online) System size dependence of pairing corre- 
lations P2d( r ) at U/t = 6 and 5 ~ 0.2. Triangles, squares and 
circles represent the 10 x 10, 8 x 8 and 6x6 lattices. The pairing 
correlations at U = are illustrated by the crosses (for 20 x 20 
lattice). The dashed line represents the asymptotic r~ 3 scaling 
for the noninteracting system. 



each other, a common typical result in these studies is 
that (A2d) in the corresponding definition by Eq.(8) has 
the amplitude around 0.1 7 ) or even larger, 8 ) which leads 
to \P 2d (r -> oo)| > 0.01. In fact, the VMC study 9 ) indi- 
cates \P2d(r —> oo ) | — 0.01 — 0.02 consistently. However, 
these are one order of magnitude larger than our upper 
bound. 

The overestimate of the order parameter by the DCA, 
CDMFT, VCT and VMC may be ascribed to the en- 
hancement of superconducting correlations found only 
in the short-ranged part in the present study. In fact, it 
is natural to expect that the enhancement in the short- 
ranged part generates overestimated pairing in the mean- 
field type treatment or in the small cluster studies in 
which long-ranged part of the fluctuations is ignored. It 
would be desired to evaluate the long-range part of corre- 
lations in extensions of the dynamical mean field theory 
to make the size scaling to the thermodynamic limit pos- 
sible and take into account long-range fluctuation effects 
in the VMC studies. 

The superconducting order parameter A2d has close 
connection to the estimate of other physical quantities 
such as the condensation energy Q, the single-particle 
gap A2d and the superconducting transition temperature 
T c . If |P2d| would be one order of magnitude reduced, 
since Q may be scaled by \P2d(r —> oo)| = A| d , the con- 
densation energv obtained as Q ~ 10" 3 t 9 ' 15 ) should be 
reduced to c± 10 _4 t. This upperbound of Q with t ~ 0.4 
eV, which is a realistic value of the cuprate superconduc- 
tors, results in Q < 0.04 meV. This is at least one or- 
der of magnitude smaller than the experimental results 
(-0.2-0.4meV). 16 18 ^ 

The single-particle gap amplitude A2d in the antin- 
odal direction has been estimated by CDMFT to be 
around O.lt. 19 ^ However, this gap is scaled linearly by 
the superconducting order parameter A2d- When the 
overestimate by the factor 3 for A2d is corrected, we 
find A2d — 0.03£ as an upperbound. If we again take 
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t ~ 0.4 eV and employ the experimentally observed gap 
amplitude A2d ~ 40-50 meV, 20 ^ the experimental value 
of — 0. It is at least factor 3 larger. 

Maier et al. claimed that from the DCA results, the 
transition to the 4 2 - y 2 -wave superconductor occurs at 
T c « 0.023t at S = 0.1, even at U/t = 4. 5 ) (Note 
that "T c " should be regarded as the critical tempera- 
ture of the Beresinski-Kosterlitz-Thouless transition if a 
strictly two-dimensional system is considered, whereas it 
is known that the additional three dimensionality would 
not change T c much.) Since T c is scaled by the gap am- 
plitude, it may again be scaled by \f\P2dix ~^ an d 
T c may be overestimated by DCA at least by a factor 3. 
Again by taking t ~ 0.4eV, the upperbound of T c may 
be 30 K, which does not offer a realistic model for the 
cuprate high-T c superconductors with T c > 100K. 

For the moment, it is not known whether this restric- 
tion could be relaxed for larger U beyond It. At least in 
the intermediate coupling region with U comparable to 
the bare bandwidth, which is realistic for the cuprate su- 
perconductors, the Hubbard model does not seem to offer 
a realistic account of the superconductivity in the right 
order of amplitude. The form factor / spatially more ex- 
tended than Eq.(9) does not change these difficulties. 

Effects of reduced quasiparticle weight on the pair- 
ing 21 ) is an unexplored and important issue to be clarified 
for theoretical prediction on the relation between P2d and 
the single-particle gap as well as T c and Q. However, it 
does not alter our conclusion that the available results on 
the Hubbard model mentioned above 5 9 ) and t-J mod- 
els 22 ' 23 ) obtained by calculating the same order param- 
eters from variational or mean-field approximations do 
not provide a clue for understanding the high T c super- 
conductivity because of the substantial overestimates of 
the pairing correlations in their approaches. Accuracy of 
numerical approaches higher than that obtained by these 
methods is required. 

We have calculated the superconducting correlations 
of the doped Hubbard model by using the PR-GBMC 
method. The present result poses constraints on the oc- 
curence of the superconductivity in the Hubbard model 
on the square lattice. Up to U/t = 7, the d-wave cor- 
relation is smaller than 10 -3 at long distances, which is 
much smaller than the recent approximate estimates by 
the variational Monte Carlo methods and extensions of 
the dynamical mean field theory. Comparing with avail- 
able theoretical works we estimated the upper bounds 
for the amplitudes of the superconducting gap, conden- 
sation energy and T c . The present estimates of the upper 
bounds are far below those of the copper oxide supercon- 
ductors. Thus our conclusion is that the simple Hubbard 
model does not offer a model for the superconductiv- 
ity in the right order of amplitude. For further progress, 
we need an accurate estimate of the correlations at long 
distances by fully taking account of fluctuations. Our re- 
sults show that the short-ranged Coulomb repulsion by 
itself does not automatically guarantee the emergence 
of the superconductivity at high critical temperatures 



even when it is close to the antiferromagnetic order in 
two dimensional systems or in the close proximity to the 
Mott insulator. This is consistent with the fact that many 
doped Mott insulators and metals near the antiferromag- 
netic quantum critical point do not automatically show 
the superconductivity even when the residual resistivity 
is very low. 25 ) We are urged to examine more detailed 
conditions beyond the simple framework of the Hubbard 
model, such as the internal structure of the Mott criti- 
cality 24 ) to understand the cuprate superconductors. 

A part of our computation has been done at the su- 
percomputer center at ISSP, Univ. of Tokyo. This work 
is partially supported by a Grant-in- Aid under the grant 
numbers 17071003 and 16076212 from MEXT, Japan. 
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